clear;
clc;close all;filebrowser;workspace;warning('off','all');
%% Set flags
massFlag = 1;
%%-----MCP Parameters--------------
TMCPacquire = 2*10^5; %set in the MCP acquisition software, in [ns]
dtMCPres = 0.025117; %TDC resolution [ns]
MCPRadius = 42;     %MCP detector radius [mm]
%%-----Parameters for data analysis--------
Xc = -10.250*0;      %[mm] center for plotting in detector frame
Yc = 8.75*0;    %[mm] center for plotting in detector frame
Angle = -30*0;   %[deg] rotation angle for plotting
tbinsize = 5; %20 [ns] ignore this if you are histcounts mass
Nbin = floor(TMCPacquire/tbinsize);
tEdge = (0:1:Nbin)*tbinsize; %[ns] for TOF binning
tList = (0.5:1:(Nbin-0.5))*tbinsize/1000; %[us]
% f(us) = 0.0214+ 2.481*sqrt(mass(amu));  %This formula from Yu's calibration
% massList = (tList/1000./15.69).^2*40-0.26;  %[amu] mass
if massFlag
    dmass = 0.25;
    massEdge = (20-dmass/2):dmass:(300+dmass/2);    %[amu]
    massList = 20:dmass:300;        %[amu]
    %     tList = 21.4 + 2481*sqrt(massList);  %[ns]
end
dXY = 0.15;        %was 0.5 [mm] spatial resolution
Xmin = -45;     %[mm] -20 plot range
Ymin = Xmin;    %[mm]
Xmax = 45;      %[mm] 20
Ymax = Xmax;    %[mm]
Xedge = Xmin:dXY:Xmax;
Yedge = Ymin:dXY:Ymax;
Xplot = (Xmin+dXY/2):dXY:(Xmax-dXY/2);
Yplot = (Ymin+dXY/2):dXY:(Ymax-dXY/2);
LabelFontSize = 12;
dmTOF = 6;      %% mass range for removing dominant peak in mass spectrum
mlist1 = [40 87 127];
species1 = {'K^+' 'Rb^+' 'KRb^+'};
dmPlot = 6;     %% mass range for plotting [amu]
mlistPlot = [127 87 40 80 174 214 167 254];     %for spatial plot
speciesPlot = {'KRb_1^+' 'Rb^+' 'K^+' 'K_2^+' 'Rb_2^+' 'KRb_2^+' 'K_2Rb^+' 'K_2Rb_2^+'};%for spatial plot
colorbg = zeros(1,3);       %Black
% colorbg = ones(1,3);        %white
%%---------Define folderpath-----------
% folderpath = 'C:\Users\KRb\lmf2txt';
% folderpath = '/Users/MGHu/Documents/KRb reaction/Data/Figure 3 data';
%% ------load .dat files--------------------------
%%%--New laser at 1kHz repetition rate with 170us ODT off before UV
filenameList = {'20190630_285nm_170us_ODT_off'};
NexpCycle = 950; % Not MCP trigger
NUV_signal = 500;   %UV pulse number for product signal per cycle
NUV_bkgd = 500;     %UV pulse number for background per cycle
t0 = 40;     %[ms]  time delay between STIRAP and 1st UV pulse

if 1
    Nfile = length(filenameList);
    realHitList = cell(1, Nfile);
    NUVpulse = zeros(1,Nfile);
    dataTOFList = cell(1, Nfile);
    dataMassList = cell(1, Nfile);
    dataXYAllList = cell(1, Nfile);
    dataXYPlotList = cell(1, Nfile);
    for i = 1:Nfile
        filename = filenameList{i};
        if (~exist([filename,'.dat'],'file'))
            error(['Please use mainMCPtxtAnalysis.m to convert /lmf2txt/',filename,...
                '.lmf.txt to ' filename, '.dat']);
        else
            realHitList{i} = dlmread([filename,'.dat'],'\t');
            %1st column ---- TOF [ns]
            %2nd column ---- X position [mm]
            %3rd column ---- Y position [mm]
            %4th column ---- MCP trigger #
            %5th column ---- Exp cycle #
        end
        NUVpulse(i) = max(realHitList{i}(:,4));
        dataTOF = zeros(NUVpulse(i),length(tList));
        dataMass = zeros(NUVpulse(i), length(massList));
        xtAll = realHitList{i}(:,1);
        XAll = realHitList{i}(:,2) - Xc;
        YAll = realHitList{i}(:,3) - Yc;
        dataXYAll = cell(1, NUVpulse(i));
        dataXYi = cell(1, NUVpulse(i));
        for j = 1:NUVpulse(i)
            xt = xtAll(realHitList{i}(:,4) == j); %[ns]
            dataTOF(j,:)= histcounts(xt, tEdge);
            xmass = 0.16248*(xt/1000+0.002891).^2; %[amu] ((xt-21.4)/2481).^2;  %[amu]
            dataMass(j,:) = histcounts(xmass, massEdge);
            %--All species at UVpulse j----
            dataX0 = XAll(realHitList{i}(:,4) == j);
            dataY0 = YAll(realHitList{i}(:,4) == j);
            dataX = dataX0*cos(-Angle*pi/180)-dataY0*sin(-Angle*pi/180);
            dataY = dataX0*sin(-Angle*pi/180)+dataY0*cos(-Angle*pi/180);
            dataXYAll{j} = histcounts2(dataY,dataX,Yedge,Xedge);
            dataXYj = zeros(length(Yplot),length(Xplot),length(mlistPlot));
            %--selected species at UVpulse j--------
            for k = 1:length(mlistPlot)
                dataRbX = dataX((xmass>=(mlistPlot(k)-dmPlot/2))&(xmass<=(mlistPlot(k)+dmPlot/2)));
                dataRbY = dataY((xmass>=(mlistPlot(k)-dmPlot/2))&(xmass<=(mlistPlot(k)+dmPlot/2)));
                dataXYj(:,:,k) = histcounts2(dataRbY,dataRbX,Yedge,Xedge);
            end
            dataXYi{j} = dataXYj;
        end
        dataTOFList{i} = dataTOF;       %[ns]
        dataXYAllList{i} = dataXYAll;      %[mm]
        dataXYPlotList{i} = dataXYi;
        dataMassList{i} = dataMass;     %[amu]
    end
    disp('UV pulse per cycle are:');
    disp(NUVpulse);
    disp('Experiment cycle is');
    disp(num2str(NexpCycle));
end

%% figures
h1 = figure('units','normalized','WindowStyle','docked','outerposition',[0 0 1 1]);
h2 = figure('units','normalized','WindowStyle','docked','outerposition',[0 0 1 1]);
h3 = figure('units','normalized','WindowStyle','docked','outerposition',[0 0 1 1]);
h4 = figure('units','normalized','WindowStyle','docked','outerposition',[0 0 1 1]);
h5 = figure('units','normalized','WindowStyle','docked','outerposition',[0 0 1 1]);
h6 = figure('units','normalized','WindowStyle','docked','outerposition',[0 0 1 1]);
h7 = figure('units','normalized','WindowStyle','docked','outerposition',[0 0 1 1]);
h8 = figure('units','normalized','WindowStyle','docked','outerposition',[0 0 1 1]);


%%----Plot total Mass spectrometer-----------------
if 1
    figure(h1);
    dataMassSum = zeros(size(dataMassList{1}(1,:)));
    for i=1:Nfile
        for j=1:500
            dataMassSum = dataMassSum+dataMassList{i}(j,:);
        end
    end
    mTOF = massList;
    countsTOF1 = dataMassSum;
    subplot(2,1,1)
    plot(massList, countsTOF1,'Color','black');
    x_product_plot_min = 20;
    x_product_plot_max = 270;
    xlim([x_product_plot_min, x_product_plot_max]);
    xlabel('m/Z');
    ylabel('Ion (count)');
    %     title(['Mass spectrum of all shots']);
    xt = get(gca, 'XTick');
    set(gca, 'FontSize', LabelFontSize);
    xt = get(gca, 'YTick');
    set(gca, 'FontSize', LabelFontSize);
    subplot(2,1,2)
    for i = 1:length(mlist1)
        countsTOF1(mTOF>=(mlist1(i)-dmTOF/2)& mTOF<=(mlist1(i)+dmTOF/2)) = 0;
    end
    plot(massList, countsTOF1,'Color','black');
    hold on
    xlim([x_product_plot_min, x_product_plot_max]);
    y_product_plot_max = 300;
    ylim([0, y_product_plot_max]);
    mlist10 = [40 87 127];
    for i = 1:length(mlist10)
        xdata = mTOF(mTOF>=(mlist10(i)-dmTOF/2)&mTOF<=(mlist10(i)+dmTOF/2));
        ydata = dataMassSum(mTOF>=(mlist10(i)-dmTOF/2)&mTOF<=(mlist10(i)+dmTOF/2));
        plot(xdata, ydata, '-', 'LineWidth', 0.5, 'Color', [0 0 0]+0.8);
        hold on;
    end
    %     dmTOF = 11;
    %     mlist = [80 174 214 254 167];
    %     species = {'K_2^+' 'Rb_2^+' 'KRb_2^+' 'K_2Rb_2^+' 'K_2Rb'};
    %     colorlist = {'red' 'cyan' 'magenta' 'blue' 'green'};
    mlist = [40 87 127 80 174];
    species = {'K^+' 'Rb^+' 'KRb^+' 'K_2^+' 'Rb_2^+'};
    colorlist = {'k' 'k' 'k' 'red' 'b'};
    for i = 1:length(mlist)
        xdata = mTOF(mTOF>=(mlist(i)-dmTOF/2)&mTOF<=(mlist(i)+dmTOF/2));
        ydata = dataMassSum(mTOF>=(mlist(i)-dmTOF/2)&mTOF<=(mlist(i)+dmTOF/2));
        plot(xdata, ydata, '-', 'LineWidth', 0.75, 'Color', colorlist{i});
        hold on;
        disp(['Numb of ', species{i},'=',num2str(sum(ydata))]);
    end
    hold off
    %     title(['UV pulse time after STIRAP: ', num2str((UVpulseNoList(k)-1)*100+15), ' ms']);
    xlabel('m/Z');
    ylabel('Ion (count)');
    xt = get(gca, 'XTick');
    set(gca, 'FontSize', LabelFontSize);
    xt = get(gca, 'YTick');
    set(gca, 'FontSize', LabelFontSize);
end

%%----Plot total Mass spectrometer with and without KRb -----------------
if 1
    figure(h2);
    dataMassSum = zeros(size(dataMassList{1}(1,:)));
    
    subplot(2,1,1)
    %%---New added-------
    dataMassSumbkgd = zeros(size(dataMassList{1}(1,:)));
    for i=1:Nfile
        for j= NUV_signal+1:NUV_signal+NUV_bkgd
            dataMassSumbkgd = dataMassSumbkgd+dataMassList{i}(j,:);
        end
    end
    countsTOFbkgd = dataMassSumbkgd;
    plot(massList, countsTOFbkgd, 'Color', 'g', 'LineWidth', 2);
    hold on
    %%%---------------
    for i=1:Nfile
        for j=1:NUV_signal
            dataMassSum = dataMassSum+dataMassList{i}(j,:);
        end
    end
    mTOF = massList;
    countsTOF1 = dataMassSum;
    
%     for i = 1:length(mlist1)
%         countsTOF1(mTOF>=(mlist1(i)-dmTOF/2)& mTOF<=(mlist1(i)+dmTOF/2)) = -1;
%     end
    plot(massList, countsTOF1,'Color', [0 0 0]+0.6*0, 'LineWidth', 0.5);
    hold on
    xlim([x_product_plot_min, x_product_plot_max]);
    ylim([0, y_product_plot_max]);
    
    %     for i = 1:length(mlist10)
    %         xdata = mTOF(mTOF>=(mlist10(i)-dmTOF/2)&mTOF<=(mlist10(i)+dmTOF/2));
    %         ydata = dataMassSum(mTOF>=(mlist10(i)-dmTOF/2)&mTOF<=(mlist10(i)+dmTOF/2));
    %         plot(xdata, ydata,'-','LineWidth',1, 'Color', 'k');
    %         hold on;
    %     end
    for i = 1:length(mlist)
        xdata = mTOF(mTOF>=(mlist(i)-dmTOF/2)&mTOF<=(mlist(i)+dmTOF/2));
        ydata = dataMassSum(mTOF>=(mlist(i)-dmTOF/2)&mTOF<=(mlist(i)+dmTOF/2));
        plot(xdata, ydata,'-','LineWidth', 1.0, 'Color', colorlist{i});
        hold on;
        disp(['Numb of ', species{i},'=',num2str(sum(ydata))]);
    end
    
    %     title(['UV pulse time after STIRAP: ', num2str((UVpulseNoList(k)-1)*100+15), ' ms']);
    xlabel('m/Z');
    ylabel('Ion (count)');
    xt = get(gca, 'XTick');
    set(gca, 'FontSize', LabelFontSize);
    xt = get(gca, 'YTick');
    set(gca, 'FontSize', LabelFontSize);
    
    hold off
    
    
    subplot(2,1,2)
    dataMassSum = zeros(size(dataMassList{1}(1,:)));
    for i=1:Nfile
        for j= NUV_signal+1:NUV_signal+NUV_bkgd
            dataMassSum = dataMassSum+dataMassList{i}(j,:);
        end
    end
    
    countsTOF1 = dataMassSum;
    
    for i = 1:length(mlist1)
        countsTOF1(mTOF>=(mlist1(i)-dmTOF/2)& mTOF<=(mlist1(i)+dmTOF/2)) = 0;
    end
    plot(massList, countsTOF1,'Color','black');
    hold on
    xlim([x_product_plot_min, x_product_plot_max]);
    ylim([0, y_product_plot_max]);
    for i = 1:length(mlist10)
        xdata = mTOF(mTOF>=(mlist10(i)-dmTOF/2)&mTOF<=(mlist10(i)+dmTOF/2));
        ydata = dataMassSum(mTOF>=(mlist10(i)-dmTOF/2)&mTOF<=(mlist10(i)+dmTOF/2));
        plot(xdata, ydata,'-','LineWidth',1, 'Color', [0 0 0]+0.8);
        hold on;
    end
    for i = 1:length(mlist)
        xdata = mTOF(mTOF>=(mlist(i)-dmTOF/2)&mTOF<=(mlist(i)+dmTOF/2));
        ydata = dataMassSum(mTOF>=(mlist(i)-dmTOF/2)&mTOF<=(mlist(i)+dmTOF/2));
        plot(xdata, ydata,'-','LineWidth', 1.0, 'Color', colorlist{i});
        
        hold on;
        disp(['Numb of ', species{i},'=',num2str(sum(ydata))]);
    end
    hold off
    %     title(['UV pulse time after STIRAP: ', num2str((UVpulseNoList(k)-1)*100+15), ' ms']);
    xlabel('m/Z');
    ylabel('Ion (count)');
    xt = get(gca, 'XTick');
    set(gca, 'FontSize', LabelFontSize);
    xt = get(gca, 'YTick');
    set(gca, 'FontSize', LabelFontSize);
end

%%-----Plot All XY distribution of selected species--------------
if 1
    figure(h3);
    for i = 1:length(mlistPlot)
        dataXYPlot = zeros(length(Yplot),length(Xplot));
        for j = 1:Nfile
            for k = 1: NUV_signal
                dataXYPlot = dataXYPlot+dataXYPlotList{j}{k}(:,:,i);
            end
        end
        subplot(3,3,i)
        maxCounts = ceil(max(max(dataXYPlot)));
        dataXYPlotCut = dataXYPlot;
        if i==2
            dataXYPlotCut(dataXYPlotCut > maxCounts/160) = maxCounts/160;
        end
        %         if i==3
        %             dataXYPlotCut(dataXYPlotCut > maxCounts/30) = maxCounts/30;
        %         end
        if maxCounts < 2
            dataXYPlotCut(1,1) = 2;
        end
        imagesc(Xplot, Yplot, dataXYPlotCut);
        if i == 1
            cmap = jet(maxCounts);
        end
        cmap(1:1,:) = colorbg;    % Make values 0 black:
        colormap(cmap);
        colorbar
        xlabel('x (mm)');
        ylabel('y (mm)');
        title(['All ',speciesPlot{i},': m=',num2str(mlistPlot(i)), '\pm', num2str(dmPlot/2), ' amu']);
        set(gca,'YDir','normal');
        xtick = get(gca, 'XTick');
        set(gca, 'FontSize', LabelFontSize);
        ytick = get(gca, 'YTick');
        set(gca, 'FontSize', LabelFontSize);
        hold on
        t = linspace(0,2*pi);
        plot(MCPRadius.*cos(t),MCPRadius.*sin(t),'w--');
    end
end

%%-----plot selected species--------
if 1
    figure(h4);
    LabelFontSize = 15;
    i = 4; %mlistPlot = [127 87 40 80 174 214 167 254];     %for spatial plot
    dataXYPlot = zeros(length(Yplot),length(Xplot));
    for j = 1:Nfile
        for k = 1: NUV_signal
            dataXYPlot = dataXYPlot+dataXYPlotList{j}{k}(:,:,i);
        end
    end
    maxCounts = ceil(max(max(dataXYPlot)));
    cmap = jet(maxCounts);
    dataXYPlotCut = dataXYPlot;
    if i == 1
        xplc = -9.975;
        yplc = 8.475;
    end
    if i==2
        dataXYPlotCut(dataXYPlotCut > maxCounts/200) = maxCounts/200;
    end
    if i==3
        dataXYPlotCut(dataXYPlotCut > maxCounts/20) = maxCounts/20;
    end
    if i==4
        xplc = -9.075;
        yplc = 7.575;
    end
    if i==5
        xplc = -10.43;
        yplc = 8.925;
    end
    if maxCounts < 2
        dataXYPlotCut(1,1) = 2;
    end
    imagesc(Xplot, Yplot, dataXYPlotCut);
    cmap(1:1,:) = colorbg;    % Make values 0 black:
    colormap(cmap);
    colorbar
    xlabel('x (mm)');
    ylabel('y (mm)');
    title([speciesPlot{i},': m=',num2str(mlistPlot(i)), '\pm', num2str(dmPlot/2), ' amu']);
    set(gca,'YDir','normal');
    xtick = get(gca, 'XTick');
    set(gca, 'FontSize', LabelFontSize);
    ytick = get(gca, 'YTick');
    set(gca, 'FontSize', LabelFontSize);
    hold on;
    t = linspace(0,2*pi);
    plot(MCPRadius.*cos(t),MCPRadius.*sin(t),'w-');
    
    axes('Position',[0.46 0.21 0.29 0.35])
    box on;
    imagesc(Xplot, Yplot, dataXYPlotCut);
    xlabel(' ');
    ylabel(' ');
    set(gca,'YDir','normal');
    dx = 8;
    xlim([xplc-dx/2,xplc+dx/2]);
    ylim([yplc-dx/2,yplc+dx/2]);
    xtick = get(gca, 'XTick');
    LabelFontSize1 = 12;
    set(gca, 'FontSize', LabelFontSize1)
    % ytick = get(gca, 'YTick');
    set(gca, 'FontSize', LabelFontSize1);
    set(gca, {'XColor', 'YColor'}, {'w', 'w'});
    cmap(1:1,:) = colorbg;    % Make values 0 black:
    colormap(cmap);
    % colorbar
    hold on
    t = linspace(0,2*pi);
    KE = 1.23984;       % [meV] 10 cm^-1
    VR = 1000;          % [V] repeller voltage
    AVMI = 46.7;        %[mm/sqrt(meV/V)]
    R10cm = AVMI*sqrt(KE/VR);
    plot(R10cm.*cos(t)+xplc, R10cm.*sin(t)+yplc, 'y--','LineWidth', 1.5);
    hold off
end

%%-----Calibrate TOF and mass relation-----------------------
if 1
    figure(h5);
    subplot(1,2,1)
    yMass = [40    87     127    87*2];
    xTOF = [15.688 23.137 27.953 32.722];
    plot(xTOF, yMass, 'o');
    MassTOF = 'A*(x+x0)^2';
    startPoints = [0.16234 0];
    fitres1 = fit(xTOF', yMass', MassTOF, 'Start', startPoints);
    hold on
    plot(fitres1, '-r');
    legend off
    hold off
    xlabel('TOF (us)');
    ylabel('Mass (amu)');
    title(['Mass(amu) =', num2str(fitres1.A), '*(TOF(us)+',num2str(fitres1.x0), ')^2']);
    subplot(1,2,2)
    plot(yMass, xTOF,  'o');
    TOFMass = 'B*sqrt(x)-y0';
    startPoints = [2.48 0];
    fitres2 = fit(yMass', xTOF', TOFMass, 'Start', startPoints);
    hold on
    plot(fitres2, '-r');
    legend off
    hold off
    xlabel('Mass (amu)');
    ylabel('TOF (us)');
    title(['TOF (us) =', num2str(fitres2.B), '*sqrt(Mass(amu))-',num2str(fitres2.y0)]);
end

%%----Plot total TOF spectrometer-----------------
if 1
    figure(h6);
    dataTOFSum = zeros(size(dataTOFList{1}(1,:)));
    
    for i=1:Nfile
        for j=1:3%NUVpulse(i)
            dataTOFSum = dataTOFSum+dataTOFList{i}(j,:);
        end
    end
    mTOF = tList;
    countsTOF1 = dataTOFSum;
    subplot(2,1,1)
    plot(tList, countsTOF1,'Color','black');
    xlim([11,35]);
    xlabel('TOF (us)');
    ylabel('Ion Count');
    title(['TOF spectrum of ions']);
    xt = get(gca, 'XTick');
    set(gca, 'FontSize', LabelFontSize);
    xt = get(gca, 'YTick');
    set(gca, 'FontSize', LabelFontSize);
    subplot(2,1,2)
    dataMassSum = zeros(size(dataMassList{1}(1,:)));
    for i=1:Nfile
        for j=1:3%NUVpulse(i)
            dataMassSum = dataMassSum+dataMassList{i}(j,:);
        end
    end
    countsTOF1 = dataMassSum;
    plot(massList, countsTOF1,'Color','black');
    xlim([20,300]);
    xlabel('Mass (amu)');
    ylabel('Ion Signal');
    %     title(['Mass spectrum of ions']);
    xt = get(gca, 'XTick');
    set(gca, 'FontSize', LabelFontSize);
    xt = get(gca, 'YTick');
    set(gca, 'FontSize', LabelFontSize);
end



%%-----Plot All XY distribution of selected species--------------
if 1
    figure(h7);
    dataTOFSum = zeros(size(dataTOFList{1}(1,:)));
    for k= 1:Nfile
        for j=1:NUV_signal
            dataTOFSum = dataTOFSum+dataTOFList{k}(j,:);
        end
    end
    j=1;
    for i = [1 5 4]    %mlistPlot = [127 87 40 80 174 214 167 254];
        dTOF =0.5; % [us]
        xTOF = 2.4807*sqrt(mlistPlot(i))-0.0016041; % [us]
        xdata = tList(tList>=(xTOF-dTOF/2)&tList<=(xTOF+dTOF/2));%[us]
        ydata = dataTOFSum(tList>=(xTOF-dTOF/2)&tList<=(xTOF+dTOF/2));
        if i == 5
            ymaxRb2 = max(ydata);
            xdataRb2 = xTOF;
        end
        if i == 1
            ymaxKRb = max(ydata);
            xdataKRb = xTOF;
        end
        if i == 4
            ymaxK2 = max(ydata);
            xdataK2 = xTOF;
        end
        subplot(2,2,j)
        plot(xdata, ydata,'o-','LineWidth', 1.0);
        if i==1
            ax = gca;
            ax.YAxis.Exponent = 3;
        end
        
        %         title([speciesPlot{i}]);
        xlim([xTOF-dTOF/2, xTOF+dTOF/2]);
        xlabel('TOF (\mus)');
        ylabel('Ion (count)');
        xticks([xTOF-dTOF/2, xTOF, xTOF+dTOF/2])
        xticklabels([xTOF-dTOF/2, xTOF, xTOF+dTOF/2])
        xt = get(gca, 'XTick');
        set(gca, 'FontSize', LabelFontSize);
        xt = get(gca, 'YTick');
        set(gca, 'FontSize', LabelFontSize);
        j = j+1;
    end
    xdata_Ring = 1:dTOF*1000; %[ns]
    ydata_Ring = zeros(size(xdata_Ring));
    ydata_Ring1 = dlmread('TOF_Lineshape_K2_Rb2.dat','\t');%simulated TOF
    ydata_Ring1 = ydata_Ring1/max(ydata_Ring1)*ymaxRb2;
    xdata_Ring = xdata_Ring + (xdataRb2-dTOF/2)*1000; %[ns]
    ydata_Ring(xdata_Ring>32615&xdata_Ring<=32615+length(ydata_Ring1)) = ydata_Ring1;
    %     xdata_Ring = ((1:length(ydata_Ring))-length(ydata_Ring)/2)/1000+32.78;  %[us]
    subplot(2,2,2)
    hold on
    plot(xdata_Ring/1000, ydata_Ring,'-r','LineWidth', 2.0);
    ylim([0, 80]);
    hold off
    
    xdata_Solid = 1:dTOF*1000; %[ns]
    ydata_Solid = zeros(size(xdata_Solid));
    ydata_Solid1 = dlmread('TOF_Lineshape_KRb.dat','\t');   %simulated TOF
    ydata_Solid1 = ydata_Solid1/max(ydata_Solid1)*ymaxKRb;
    xdata_Solid = xdata_Solid + (xdataKRb-dTOF/2)*1000; %[ns]
    ydata_Solid(xdata_Solid>27877&xdata_Solid<=27877+length(ydata_Solid1)) = ydata_Solid1;
    subplot(2,2,1)
    hold on
    plot(xdata_Solid/1000, ydata_Solid,'-r','LineWidth', 2.0);
    hold off
    subplot(2,2,3)
    hold on
    xdata_Ring = (0:sqrt(40/87):dTOF*1000); %[ns]
    ydata_Ring = zeros(size(xdata_Ring));
    ydata_Ring1 = dlmread('TOF_Lineshape_K2_Rb2.dat','\t');%simulated TOF
    ydata_Ring1 = ydata_Ring1/max(ydata_Ring1)*ymaxK2;
    xdata_Ring = xdata_Ring + (xdataK2-dTOF/2)*1000; %[ns]
    ydata_Ring(262:261+length(ydata_Ring1)) = ydata_Ring1;
    plot(xdata_Ring/1000, ydata_Ring,'-r','LineWidth', 2.0);
    ylim([0, 30]);
end

%%---- fit K2 and Rb2 VMI to 2D gaussian-------
if 1
    figure(h8);
    LabelFontSize = 12;
    iplot = 1;
    for i = 4:5    %mlistPlot = [127 87 40 80 174 214 167 254];     %for spatial plot   
        dataXYPlot = zeros(length(Yplot),length(Xplot));
        for j = 1:Nfile
            for k = 1: NUV_signal
                dataXYPlot = dataXYPlot+dataXYPlotList{j}{k}(:,:,i);
            end
        end
        maxCounts = ceil(max(max(dataXYPlot)));
        cmap = jet(maxCounts);
        dataXYPlotCut = dataXYPlot;
        if i == 1
            xplc = -9.975;
            yplc = 8.475;
        end
        if i == 4
            xplc = -9.075;
            yplc = 7.575;
        end
        if i == 5
            xplc = -10.43;
            yplc = 8.925;
        end
        if maxCounts < 2
            dataXYPlotCut(1,1) = 2;
        end
        subplot(2, 3, iplot);
        iplot = iplot + 1;
        imagesc(Xplot, Yplot, dataXYPlotCut);
        cmap(1:1,:) = colorbg;    % Make values 0 black:
        colormap(cmap);
        colorbar
        xlabel('x (mm)');
        ylabel('y (mm)');
        title([speciesPlot{i},' VMI data']);
        set(gca,'YDir','normal');
        xtick = get(gca, 'XTick');
        set(gca, 'FontSize', LabelFontSize);
        ytick = get(gca, 'YTick');
        set(gca, 'FontSize', LabelFontSize);
        dx = 8;     %[mm]
        xlim([xplc-dx/2,xplc+dx/2]);
        ylim([yplc-dx/2,yplc+dx/2]);
        hold on
        t = linspace(0,2*pi);
        KE = 10.4;       % [cm^-1]
        VR = 1000;          % [V] repeller voltage
        AVMI = 16.44;        %[mm/sqrt(cm-1/V)]
        R10cm = AVMI*sqrt(KE/VR);
        plot(R10cm.*cos(t)+xplc, R10cm.*sin(t)+yplc, 'y--','LineWidth', 1.5);
        hold off
        
        subplot(2, 3, iplot);
        iplot = iplot + 1;
        ROIxc = xplc;   %[mm]
        ROIyc = yplc;   %[mm]
        XplotCut = Xplot(Xplot>=ROIxc-dx/2&Xplot<=ROIxc+dx/2);
        YplotCut = Yplot(Yplot>=ROIyc-dx/2&Yplot<=ROIyc+dx/2);
        ROIxPX = find(Xplot>=ROIxc-dx/2&Xplot<=ROIxc+dx/2);
        ROIyPX = find(Yplot>=ROIyc-dx/2&Yplot<=ROIyc+dx/2);
        dataXYPlotCut = dataXYPlot(ROIyPX, ROIxPX);
        maxCounts = ceil(max(max(dataXYPlotCut)));
        %%--fit to 2D gaussian-----
        offset=0;
        fitpar = [maxCounts ROIxc ROIyc 0.3 0.3]; %(Amp, xc, yc, wx, wy, Amp2, wx2, wy2)
        params = offset;
        x = Xplot(ROIxPX);  %[mm]
        y = Yplot(ROIyPX);  %[mm]
        [X,Y] = meshgrid(x,y);
        Fitdata = dataXYPlotCut;
        %set options for fminsearch
        options = optimset('Display','off','TolFun',1e-6,'TolX',1e-8,'MaxFun',4000, 'MaxIter', 2000);
        [fit0,~,~,~] = fminsearch(@gaussianoffset,fitpar,options,params,X,Y,Fitdata);
        A = fit0(1);
        xc = fit0(2);
        yc = fit0(3);
        wx = fit0(4);
        wy = fit0(5);
%         A2 = fit0(6);
%         wx2 = fit0(7);
%         wy2 = fit0(8);
        area1=(X-xc).^2./(2*wx.^2)+(Y-yc).^2./(2*wy.^2);
%         area2=(X-xc).^2./(2*wx2.^2)+(Y-yc).^2./(2*wy2.^2);
        FitResult=offset + A.*exp(-area1);% + A2.*exp(-area2);
        imagesc(XplotCut, YplotCut, FitResult);
        xlabel('x (mm)');
        ylabel('y (mm)');
        title([speciesPlot{i}, ' fit: xc=', num2str(xc), ' mm, yc=', num2str(yc), ' mm']);
        set(gca,'YDir','normal');
        xtick = get(gca, 'XTick');
        set(gca, 'FontSize', LabelFontSize);
        ytick = get(gca, 'YTick');
        set(gca, 'FontSize', LabelFontSize);
        colorbar
        %%--plot residue---
        subplot(2, 3, iplot);
        iplot = iplot + 1;
        FitResidue = FitResult - dataXYPlotCut;
        imagesc(XplotCut, YplotCut, FitResidue);
        xlabel('x (mm)');
        ylabel('y (mm)');
        title([' wx=', num2str(wx), ' mm, wy=', num2str(wy), ' mm']);
        set(gca,'YDir','normal');
        xtick = get(gca, 'XTick');
        set(gca, 'FontSize', LabelFontSize);
        ytick = get(gca, 'YTick');
        set(gca, 'FontSize', LabelFontSize);
        colorbar
    end
end



figure(h2);
% savefig(h, ['plots/summaryplot_K2Rb2_spectroscopy.fig']);
% saveas(h7, ['plots/TOFlineshape.png']);
% saveas(h4, ['plots/VMI_285nm_KRb.png']);
% saveas(h2, ['plots/MassSpec_285nm.png']);